Explicit wave-averaged primitive equations 
using a Generalized Lagrangian Mean 



Fabrice Ardhuin^, Nicolas Rascle^'^, K. A. Belibassakis 

^Centre Militaire d'Oceanographie, Service Hydrographique et Oceanographique de 

la Marine, 29609 Brest, France 

^Laboratoire de Physique des Oceans, Universite de Bretagne Occidentale, 29000 

Brest, France 

'^School of Naval Architecture and Marine Engineering, National Technical 
University of Athens, Athens, Greece 



Abstract 

The generalized Langrangian mean theory provides exact equations for general 
wave-turbulence-mean flow interactions in three dimensions. For practical applica- 
tions, these equations must be closed by specifying the wave forcing terms. Here an 
approximate closure is obtained under the hypotheses of small surface slope, weak 
horizontal gradients of the water depth and mean current, and weak curvature of the 
mean current profile. These assumptions yield analytical expressions for the mean 
momentum and pressure forcing terms that can be expressed in terms of the wave 
spectrum. A vertical change of coordinate is then applied to obtain glm2z-KKHS 
equations (55) and (57) with non-divergent mass transport in cartesian coordinates. 
To lowest order, agreement is found with Eulerian-mean theories, and the present 
approximation provides an explicit extension of known wave-averaged equations to 
short-scale variations of the wave field, and vertically varying currents only limited 
to weak or localized profile curvatures. Further, the underlying exact equations pro- 
vide a natural framework for extensions to finite wave amplitudes and any realistic 
situation. The accuracy of the approximations is discussed using comparisons with 
exact numerical solutions for linear waves over arbitrary bottom slopes, for which 
the equations are still exact when properly accounting for partial standing waves. 
For finite amplitude waves it is found that the approximate solutions are proba- 
bly accurate for ocean mixed layer modelling and shoaling waves, provided that an 
adequate turbulent closure is designed. However, for surf zone applications the ap- 
proximations are expected to give only qualitative results due to the large infiuence 
of wave nonlinearity on the vertical profiles of wave forcing terms. 

Key words: radiation stresses. Generalized Lagrangian Mean, wave-current 
coupling, drift, surface waves, three dimensions 
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1 Introduction 



Prom wave- induced mixing and enhanced air-sea interactions in deep water, to 
wave-induced currents and sea level changes on beaches, the effects of waves on 
ocean currents and turbulence are well documented (e.g. Battjes 1988, Terray 
et al. 1996). The refraction of waves over horizontally varying currents is also 
well known, and the modifications of waves by vertical current shears have 
been the topic of a number of theoretical and laboratory investigations (e.g. 
Biesel 1950, Peregrine 1976, Kirby and Chen 1989, Swan et al. 2001), and 
field observations (e.g. Ivonin et al. 2004). In spite of this knowledge and the 
importance of the topic for engineering and scientific applications, ranging 
from navigation safety to search and rescue, beach erosion, and de-biasing 
of remote sensing measurements, there is no well established and generally 
practical numerical model for wave-current interactions in three dimensions. 

Indeed the problem is made difficult by the difference in time scales between 
gravity waves and other motions. When motions on the scale of the wave 
period can be resolved, Boussinesq approximation of nearshore flows has pro- 
vided remarkable numerical solutions of wave-current interaction processes 
(e.g. Chen et al. 2003, Terrile ct al. 2006). However, such an approach still 
misses some of the important dynamical effects as it cannot represent real 
vertical current shears and their mixing effects (Putrevu and Svendsen 1999). 
This shortcoming has been partly corrected in quasi-three dimensional models 
(e.g. Haas et al. 2003), or multi-layer Boussinesq models (e.g. Lynnett and Liu 
2005). 

The alternative is of course to use fully three dimensional (3D) models, based 
on the primitive equations. These models are extensively used for investigating 
the global, regional or coastal ocean circulation (e.g. Bleck 2002, Shchepetkin 
and McWilliams 2003). An average over the wave phase or period is most 
useful due to practical constraints on the computational resources, allowing 
larger time steps and avoiding non-hydrostatic mean ffows. Wave-averaging 
also allows an easier interpretation of the model result. A summary of wave- 
averaged models in 2 or 3 dimensions is provided in table 1. 



1.1 Air-water separation 

In 3D, problems arise due to the presence of both air and water in the region 

between wave crests and troughs. Various approaches to the phase or time 
averaging of ffow properties are illustrated in figure 1 (see also Ardhuin et 
al. 2007b, hereinafter ARB2007). For small amplitude waves, one may simply 
take a Taylor expansion of mean flow properties (e.g. McWilliams et al. 2004, 
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Table 1 



Essential attributes of some general wave-current coupling theories. See list of sym- 
bols for details (table 2 at the end of the paper). Although Mellor (2003) derived 

his wave-averaged equations with spatially varying wave amplitudes, his use of flat- 
bottom Airy wave kinematics is inconsistent with the presence of bottom slopes 
(see ARB07). MRL04 stands for McWilliams et al. (2004) and NA2007 stands for 
Newberger and Allen (2007). 

hereinafter MRL04). Using a decomposition of the non-linear advection term 
in the equations of motion u- Vu = Vm^ -|- u x Vu, McWilliams et al. (2004, 
see also Lane et al. 2007) obtained a relatively simple set of equation for 
conservative wave motion over sheared currents, for a given choice of small 
parameters. These parameters include the surface slope Si = /cqOo and the ratio 
of the wavelength and scale of evolution of the wave amplitucde. Further, these 
equations were derived with a scaling corresponding to a non-dimensional 
depth koho of order 1, with ko, gq and ho typical values of the wavenumber, 
wave amplitude and water depth, respectively. These authors also assumed 
that the current velocity was of the same order as the wave orbital velocity, 
both weaker than the phase speed by a factor ei. That latter assumption may 
generally be relaxed since the equations of motion are invariant by a change of 
reference frame, so that only the current vertical shear may need to be small 
compared to the wave radian frequency, provided that the current, water depth 
and wave amplitudes are slowly varying horizontally. 

For waves of finite amplitude, a proper separation of air and water in the 
averaged equations of motion requires a change of coordinates that maps the 
moving free surface to a level that is fixed, or at least slowly varying. This 
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is usual practice in air-sea interaction studies, and it has provided approxi- 
mate solutions to problems such as wind-wave generation or wave-turbulence 
interactions (e.g. Jenkins 1986, Teixeira and Belcher 2002) but it brings some 
complications. The most simple change of coordinate was recently proposed 
by Mellor (2003), but it appears to be impractical in the presence of a bot- 
tom slope because its accurate implementation requires the wave kinematics 
to first order in the wave slope (Ardhuin et al., 2007b, hereinafter ARB07). 

1.2 Separation of wave and current momentum fluxes 

Another approach is to use one of the two sets of exact averaged equations 
derived by Andrews and Mclntyre (1978a). Groeneweg (1999) successfully 
used the second set, the alternative Generalized Lagragian Mean equations 
(aGLM), approximated to second order in wave slope, for the investigation 
of current profile modifications induced by waves (see also Groeneweg and 
Klopman 1998, Groeneweg and Battjes 2003). This work was also loosely 
adapted for engineering use in the numerical model Delft3D (Walstra et al. 
2001). 

However, aGLM equations describe the evolution of the total flow momentum, 
which includes the wave pseudo- momentum per unit mass P. That vector 
quantity is generally close to the Lagrangian Stokes drift u"^ (see below), and 
it is not mixed by turbulenc^, unlike the mean flow momentum. Further, 
P is carried by the wave field at the group velocity, which is typically one 
order of magnitude faster than the drift velocity. Thus bundling P with the 
rest of the momentum may lead to large errors with the turbulence closure. 
Other practical problems arise due to the strong surface shear of P and u'^ 
(e.g. Rascle et al. 2006) whereas the quasi-Eulerian current is relatively uni- 
form in deep water (e.g. Santala and Terray 1992). Thus solving for the total 
momentum (including P) requires a high resoltion near the surface. Finally, 
a consistent expression of the aGLM equations with a sloping bottom and 
wave field gradients is difficult due to the divergence of vertical fluxes of mo- 
mentum (vertical radiation stresses) that must be expressed to first order in 
all the small parameters that represent the slow wave field evolution (bottom 
slope, wave energy gradients, current shears...). This same problem arises with 
Mellor's (2003) equations and is discussed in ARB07. 

The first set of GLM equations describes the evolution of the quasi-Eulerian 
current only, and, just like the decomposition of u- Vm used by MRL04, it does 
not require the evaluation of these vertical radiation stresses. These equations 
were used by Leibovich (1980) to derive the Craik-Leibovich equations that 

^ The Stokes drift is a residual velocity over the wave cycle, its mixing is not possible 
without a profound modification of the wave kinematics. 
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Fig. 1. Averaging procedures (left) and examples of resulting velocity profiles 
(right) in the case of (a) Eulerian averages (e.g. Rivero and Sanchez- Arcilla 1995, 
McWilliams et al. 2004), (b) the Generalized Lagrangian Mean (Andrews and Mcln- 
tyre 1978a), and (c) sigma transform (Mellor 2003, AJB07). The thick black bars 
connect the fixed points x where the average field is evaluated, to the displaced 
points X + ^ where the instantaneous field is evaluated. For averages in moving co- 
ordinates the points x + ^ at a given vertical level ^ are along the gray lines. The 
drift velocity is the sum of the (quasi-Eulerian) current and the wave-induced mass 
transport. In the present illustration an Airy wave of amplitude 3 m and wave- 
length 100 m in 30 m depth, is superimposed on a hypothetical current of velocity 
u{z) = —0.5 — O.Olz m/s for all z < C(x)- The current profile is not represented 
in (c) since it is not directly given in Mellor's theory, although it can obviously be 
obtained by taking the difference of the other two profiles. 



is the basis of theories for Langmuir circulations. However, in that work he 
did not attempt an explicit integration of the GLM set, and thus did not 
express the wave forcing terms from wave amplitudes or spectra. The general 
mathematical structure of the GLM equations and their conservatin properties 
are also well detailed in Holm (2002) and references therein. 



5 



Further, the GLM flow is generally divergent as the averaging operator intro- 
duces an implicit change of the vertical coordinate. This question has been 
largely overlooked by previous users of GLM theory (Leibovich 1980, Groe- 
neweg 1999). Further, in order to be implemented in a numerical model, the 
wave-induced forcing terms must be made explicit using approximate solu- 
tions for wave-induced motions and pressure. We will assume that the slowly 
varying spectrum is known, typically provided by a wave model. Given the 
degree of accuracy attained by modelled wave spectra in a wide variety of 
conditions this is generally appropriate (e.g. Herbers et al. 2000, Ardhuin et 
al. 2003, 2007, Magne et al. 2007). We note in passing that no exphcit and 
theoretically satisfying theory is available for the transport of the wave action 
spectrum over vertically and horizontally sheared currents. Indeed, the exact 
theory of Andrews and Mclntyre (1978b) is implicit and would require an ex- 
plicit approximation of the wave action from know wave kinematics, similar 
to the approximation of the wave pseudo-momentum performed here. 

The goal of the present paper is to provide a practical and accurate method for 
wave-current coupling that is general enough for applications ranging from the 
ocean mixed layer to, possibly, the surf zone. GLM equations, for the reasons 
listed above, are a good candidate for this application. Although not as simple 
as an Eulerian average, the GLM operator is capable of properly separating 
air and water in the crest to trough region, leading to physically understand- 
able definitions of mean properties on either side of the air-sea interface. The 
practical use of GLM requires some approximations and transformations. We 
provide in section 2 a derivation of explicit and approximate gl'm2z-RANS 
equations. Given the large literature on the subject, we explore in section 3 
the relationships between GLM, aGLM and other forms of wave- averaged 3D 
and depth-integrated 2D equations. A preliminary analysis of the expected 
errors due to the approximations are provided in section 4, and conclusions 
follow in section 5. Full numerical solutions using the glm2z-RANS equations 
will be reported elsewhere, in particular in the doctorate thesis of Nicolas 
Rascle. 



2 glm2-RANS equations 

2.1 Generalities on GLM and linear wave kinematics 



We first define the Eulerian average (j) (x, t) of (j) (x, t), where the average may 
be an average over phase, realizations, time t or space. We now take this 
average at displaced positions x+^, with ^ = (^i, ^3) a displacement vector, 
and we defining the velocity v at which the mean position is displaced when 
the actual position moves at the fiuid velocity u(x + ^). One obtains the 
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corresponding GLM of 



0(x,t)'' = 0(x+e,t) (1) 

by choosing the displacement field ^ so that 

• the mapping x — > x + ^ is invertible 

. eM) = ^ 

• V (x, t) = V (x, t), which gives v = u(x, t) . 

Such a mapping is illustrated in figure l.c for linear waves. Lagrangian per- 
turbations are logically defined as the field minus its average, i.e., 

0(x, tf = 0(x + ^,t)- 0(x, tf = 0(x + t) - 0(x + ^,t). (2) 

Here we shall take our Eulerian average to be a phase averagj^. Given any 
Eulerian fiow field u(x, t), one may define a first displacement by 

t+At 

^'(x, t, At) = 1 u(x + r(x, t, t' - t),t')dt'. (3) 

t 



The mean drift velocity is defined as v(x, t) = limAt^o ^'(x, t, At)/(At). The 
GLM displacement field is then given by ^ = vt— ^' — vt. This construction 
of V and ^ guarantees that the required properties are obtained, provided that 
the limit At commutes with the averaging operator. For periodic motions 
one may also take v = (^'(t + T^) — ^'(t))/ (T^), with the Lagrangian wave 
period (the time taken by a water particle to return to the same wave phase). 
This definition will be used for Miche waves in section 4.2. 

Clearly GLM differs from the Eulerian mean. The difference between the two 
is given by the Stokes correction (Andrews et Mclntyre 1978a). Below the 
wave troughs, the Stokes correction for the velocity is the Stokes drift, by 
definition, 

= - u. (4) 



^ For uncorrelated wave components the phase average is obtained by the sum of 
the phase averages of each component. In the presence of phase correlations, such 
as in the case of partially standing waves or nonlinear phase couplings, the sum has 
to be averaged in a coherent manner. 



7 



More generally, for a continuously differentiablc field the Stokes correction 
is given by (Andrews and Mclntyre 1978a, equation 2.27), 



with an implicit summation over repeated indices. 

The GLM average commutes with the Lagrangian derivative, thus the GLM 
velocity is the average drift velocity of water particles. One should however 
be careful that the GLM average does not commute with most differential op- 
erators, for example the curl operator. Indeed the GLM velocity of irrotational 
waves is rotational, which is clearly apparent in the vertical shear of the Stokes 
drift (see also Ardhuin and Jenkins 2006 for a calculation of the lowest order 

mean shears dua/dz and du^/dx ). 

One of the interesting aspects of GLM theory is that it clearly separates the 
wave pseudo-momentum P from the quasi-Eulerian mean momentum u = 
u-^ — P. This is a key aspect for numerical modelling since P is transported by 
the wave field at the group velocity, of the order of 5 m s~^ in deep water, while 
u is transported at the much slower velocity u^. P is defined by (Andrews 
and Mclntyre 1978a, eq. 3.1), 



where CijkAjBk is the i-component of the vector product A x B, and /fe/2 is the 
/c-component of the rotation vector of the reference frame. In the applications 
considered here the effect of rotation can be neglected in (6) due to the much 
larger rotation period of the Earth compared to the wave period. We will thus 
take 



For practical use, the GLM equations have to be closed by specifying the 
wave-induced forcing terms. In order to give explicit approximations for the 
wave-induced effects, we will approximate the wave motion as a sum of linear 
wave modes, each with a local wave phase ip giving the local wave number 
k = {ki, — VV', and radian frequency cu — —dip/dt, and an intrinsic linear 
wave radian frequency a — [gktanh.{kD)]^^^ — oj — k-L/^, where XJ a is the 
phase advection velocity, D is the local mean water depth, and g the accelera- 
tion due to gravity and Earth rotation. Defining h{xi,X2) as the local depth of 
the bottom and Ci^i^ ^2, t) as the free surface elevation, one has D = ( + h. We 
assume that the wave slope £1 = max (| VC|) is small compared to unity (this 




(6) 




(7) 
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will be our first hypothesis HI), with V denoting the horizontal gradient op- 
erator. We also restrict our investigations to cases for which the Ursell number 
is small Ur = {a/D)/ {kDy < 1 (this is hypothesis H2). We further restrict 
our derivations to first order in the slow spatial scale £2- That small parame- 
ter may be defined as the maximum of the slow spatial scales \(da/dx)/{ka)\, 
\{du/dx)/{a)\, |(9D/(9a;)|, and time scales \{da/dt)/{aa)\, \k{du/dt) /{a)"^, a,n.d 
\{dD/dt)k/a\ (hypothesis H3). It will also appear that the current profile may 
cause some difficulties. Since we have already assumed a small wave steepness 
we may use Kirby and Chen's (1989) results, giving the dispersion relation 



, 2kcosh\2k(z + h)] ^ , 



where a is a dummy index representing any horizontal component 1 or 2, and 
the summation is implicit over repeated indices. The index 3 will represent the 
vertical components positive upwards, along the direction z — Xs.ln particular 
we shall assume that their correction to the lowest order stream function 
(their eq. 23) is relatively small, which may be obtained by requiring that the 
curvature of the current is weak or concentrated in a thin boundary layer, i.e. 
£3 -C 1 (hypothesis H4) with 



£3 = 



1 



u;smh.{kD) 



sinh [2k{z + h)] dz. 



(9) 



For simplicity we will further require that [d^Ua/ dz^ / {a)] < £3 (hypothe- 
sis H5), which may be more restrictive than H4. Finally, we will neglect the 
vertical velocity w in the vertical momentum equation for the mean flow mo- 
mentum (i.e. we assume the mean flow to be hydrostatic, this is our hypothesis 
H6). 

In the following we take e = max£j, 1 ^ i ^ 3. The wave-induced pressure and 
velocity are given by 



p = Py,ga[Fcccosip + 0{s)] (10) 
k 

aa^ [Fes cos ijj + 0(s)] (11) 
k 

U3^aa[Fsssmil; + 0{e)], (12) 

where a is the local wave amplitude, is the water density, taken constant 
in the present paper. We have used the short-hand notations Fee = cosh{kz + 
kD)/cosh{kD), Fes = cosh{kz+kD)/ smh{kD), and Fss = smh{kz+kD)/ smh{kD). 
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From now on, only the lowest order approximations will be given unless explic- 
itly stated otherwise. In order to estimate quantities at displaced positions, 
the zero-mean displacement field is given by 



Thanks to the definition of u^, we also have 

dt dx, ~ dt '^dx' 



I _Ldii dii _^ dii 

= ^ + % ^ - ^ + ^ ' (14) 



in which the vertical velocity has been neglected. The greek indices a and (3 
stand for horizontal components only. 

To lowest order in the wave amplitude, the displacements ii and Lagrangian 
velocity perturbations u\ are obtained from (13) and (14), 



$,3 = am [Fss cos^] 



u'^^Ua + 6^ + ^pp^ + O (aka'') cos 2^^ + O ( a" 



dz 



dx 



dz^ 



ia = -am 



ka „ m dua „ 

-rFcs H 'E~^ss 

k a dz 



^ , a dua \ , ^ 
+0 { 1 cos^ O 



sin -|- O 



a dz"^ 



sin 2ip 



a dx 



dz^ 



(15) 
(16) 

(17) 
(18) 



(19) 



The shear correction parameter m, arising from the time-integration of (14), 
is given by 



m(x, k, z, t) 



uj — k'U (x, z^ t) 



(20) 



Based on (8) m differs from 1 by a quantity of order a ^dujdz. 

Using our assumption (H5) the last term in eq. (19) may be neglected. The 
last two term in eq. (17) have been neglected because they will give negligible 
0(£^) terms in P, ^ or other wave-related quantities, when multiplied by 
other zero-mean wave quantities. 
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Using the approximate wave-induced motions, one may estimate the Stokes 
drift 



2 



ma 



4sinh'(fcD) 

sinh^(A;^ + kh) 



2crkcosh(2A;2; + 2kh) + kmsinh(2A;2; + 2kh) 



k dn 

k dz 



+ 



the horizontal wave pseudo-momentum 



(21) 



ma 



4sinh^(A;D) 

k / du^ 

+2m^-^smh'^{kz + kh) ( — 



and the GLM position of the free surface 



2akc,cosh{2kz + 2kh) + 2kc,msmh{2kz + 2kh)^- — 



- ma 



(22) 



k mk du 

tanh/cD a dz 



• (23) 



Thus the GLM of vertical positions in the water is generally larger than the 
Eulerian mean of the position of the same particles (see also Mclntyre 1988). 
This is easily understood, given that there are more particles under the crests 
than under the troughs (figure l.c). As a result, the original GLM equations 
are divergent (V-u^ ^ 0) and require a coordinate transformation to yield a 
non-divergent velocity field. That transformation is small, leading to a relative 
correction of order e\. That transformed set of equation is a modified primitive 
equation that may be implemented in existing ocean circulation models. 



The horizontal component of the wave pseudo-momentum differs from the 
Stokes drift u^^ due to the current vertical shear. Therefore the quasi- Eulerian 
mean velocity Ua = — Pa also differs from the Eulerian mean velocity 



(24) 



The vertical wave pseudo-momentum P3 = is, at most, of order ae^ /k. Al- 
though it may be neglected in the momentum equation, it plays an important 
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role in the mass conservation equation, and will thus be estimated from Pq. 
In particular, for m = 1 and in the limit of small surface slopes, it is straight- 
forward using (7) to prove that P is non-divergent and such that P-n = at 
z — —h, with n the normal to the bottom. This gives. 



Although this equality is not obvious for m 7^ 1 and nonlinear waves, correc- 
tions to (25) are expected to be only of higher order ,in particular once P is 
transformed to z coordinates. Indeed, in the absence of a mean flow P = 
and it is non-divergent (see section 2.1.1). 

2.1.1 glm2-RANS equations 

The velocity field is assumed to have a unique decomposition in mean, wave 
and turbulent components u = u -|- u -|- u', with (u') = 0, the average over 
the flow realizations for prescribed wave phases. The turbulence will be as- 
sumed weak enough so that its effect on the sea surface position is negligible. 
We note X the divergence of the Reynolds stresses, i.e. Xj = d (u'^u'^^ /dxj, 
and we apply the GLM average to the Reynolds-Average Navier-Stokes equa- 
tions (RANS). We shall now seek an approximation to the GLM momentum 
equations by retaining all terms of order p^ge^ and larger in the horizontal 
momentum equation, and all terms of order p^g^^ in the vertical momentum 
equation. The resulting equations, that may be called the " 5'/m2-RANS" equa- 
tions, are thus more limited in terms of wave nonlinearity than the Eulerian 
mean equations of MRL04. At the same time, random waves are considered 
here and that the mean current may be larger than the wave orbital velocity. 
Indeed we make no hypothesis on the current magnitude, but only on the 
horizontal current gradients and on the curvature of the current profile. The 
present derivation differs from that of Groeneweg (1999) by the fact that we 
use the GLM instead of the aGLM equations (see table 1). The name for these 
equations is loosely borrowed from Holm (2002) who instead derived an ap- 
proximate Lagrangian to obtain the momentum equation, and did not include 
turbulence. 

In order to simphfy our calculations we shall use the form of the GLM equa- 
tions given by Dingemans (1997, eq. 2.596) with p^ constant, which, among 
other things, removes terms related to the fluid thermodynamics. The evolu- 
tion equation for the quasi- Eulerian velocity u is. 




(25) 




(26) 
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where the Lagrangian derivative is a derivative following the fluid at the 
Lagrangian mean velocity U'^, p is the full dynamic pressure, 5 is Kronecker's 
symbol, and the viscous and/or turbulent force X is defined by 



These exact equations will now be approximated using (10)-(16). We first 
evaluate the wave forcing terms in (26) using monochromatic waves, with a 
surface elevation variance E = a^/2. The result for random waves follows by 
summation over the spectrum and replacing E with the spectral density E(k). 

We first consider the vertical momentum balance, giving the pressure field. It 
should be noted that the Lagrangian mean Bernoulli head term ■u^m^/2 differs 
from its Eulerian counterpart u'jU'j/2 by a term K2, which arises from the 
correlation of the mean current perturbation at the displaced position x + ^, 
with the wave-induced velocity, i.e. the second term in (17). Eqs. (10)-(16) 
give 

I (44) = ^ [FccFcs + FscFss] + i^2, (28) 



with 



K2 = Ma6 



dz 



2 



(9u 



dz 



E-k-—mFcsFss + — 



du 



dz 



ss- 



(29) 



The vertical momentum equation (26) for id = €,3 is. 



did div „ dm „ , 



did 



dt 



dz 



dz 



dx 



1 dp^ 
Pm dz 
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d_ 

dz 



(uaUa + iv^) /2 + K2 



+ Pp— {Up + Pp) + P,— {% + P3) (30) 



For small bottom slopes we may neglect the last term, but we rewrite it in 
order to compare with other sets of equations. Now using the lowest order 
wave solution (11)-(16), eq. (30) transforms to 



— -^P +Pwgz-Pw^{Fcs + Fss)-PwK: 

pw OZ |_ / 

^ ^ dm ^ d ,^ „ 



dm ^dm 

dt dz 



{up + Pp)t^ + Pa— 



dx 



■/3 



+ Pp) + {W + P3) . (31) 
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We add to both sides the depth-uniform term —a'^E{FQQ — Fgg) /2, and in- 
tegrate over z to obtain 



^ = [{z - zs) - kEFccFcs] + K, + K^- . f^f (32) 



where the hydrostatic hypothesis (H6, see above) has be made for the mean 
flow. The depth-integrated vertical component of the vortex-hke force Ki is 
defined by 

Ki^-j Pp-Q^ {up + Pp) dz' + I (Pfi) dz', (33) 



where eq. (25) has been used. The integration constant Zs is given by the 
surface boundary condition 

MO'' = -Pw9 (f -Zs- kEFccFcs - K^{t)/g) = Pa- (34) 



Using (23) we find that Zs = C, + Pa/{Pwg) — K2{{C)^)/g and (32) becomes 



^ = ^ + gkEFccFcs + K, + K2- ^t), (35) 

Pw Pw 



with the hydrostatic pressure defined equal to the mean atmospheric pres- 
sure at the mean sea surface, — Pwg{C — z) +Pa- 

Below the wave troughs the Stokes correction for the pressure (5) gives the 
Eulerian-mean pressure 

(k dvi \ 
Fes Fee + FssFsc + '^"q^^^^s^^(^'<^' 1 ■ (2^) 



Thus equation (32) gives the following relationship, valid to order ef below 
the wave troughs, between the Eulerian-mean pressure p and p^. 



P^p" - p^gkEFssFsc + Pn, \Ki - K^iC"") + | 

+p^gk{l - m)EFccFcs- (37) 



dvL 




dz 


m'FlA 
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For a spectrum of random waves, the modified pressure term that enters the 
horizontal momentum equation may be written as 

- _ PwU^U^j duf _ CJ I ^ Qshear (oo\ 

p = p ^~dz' ~ Pw'^ ' 

with the depth-uniform wave-induced kinematic pressure term 

S'=sf *fl^' dk (39) 
smh2kD ^ ^ 

k 

and a shear-induced pressure term, due to the integral of the vertical compo- 
nent of the vortex force Ki, and K2{C ), 



qshear ^ _ y ^^^^ ^^^^ 

k 

t 

-II 



dk 



k 2 



^ ^ dz' 



d^'dk. 



(40) 



Now considering the horizontal momentum equations, we rewrite (26) for the 
horizontal velocity. 



■A (S' + + P^l^ - + ^a, (41) 



Grouping all Pp terms, as in Garrett (1976 eq. 3.10 and 3.11), leads to an 
expression with the 'vortex force' ecx^AfsOJ^Pp. This force is the vector product of 
the wave pseudo-momentum P and mean flow vertical vorticity uj^. Equation 
(41) transforms to 



The vortex force is a momentum flux divergence that compensates for the 
change in wave momentum flux due to wave refraction over varying currents. 
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and includes the flux of momentum resulting from u momentum advected by 
the wave motion (Garrett 1976). 



The turbulent closure is the topic of ongoing research and will not be explicitly 
detailed here. We only note that it differs in principle from the closure of the 
aGLM equations of Groeneweg (1999), which could be extended to include 
the second term in eq. (27). A proper closure involves a full discussion of the 
distortion of turbulence by the waves when the turbulent mixing time scale 
is larger than the wave period (e.g. Walmsley and Taylor 1996, Janssen 2004, 
Teixeira and Belcher 2002). One should consider with caution the rather bold 
but practical assumptions of Groeneweg (1999) who used a standard turbu- 
lence closure to define the viscosity that acts upon the wave-induced velocities, 
or the assumption of Huang and Mei (2003) who assumed that the eddy vis- 
cosity instantaneously adjusts to the passage of waves. These effects may have 
consequences on the magnitude of wave attenuation through its interaction 
with turbulence, and the resulting vertical profile of X^- Here we only note 
that any momentum lost by the wave field should be gained by either the 
atmosphere, the bottom or the mean flow. Thus a possible parameterization 
for the diabatic source of momentum is 



dRgp , 9 I dUg ^ T^turb T^bfric rAo\ 



with Rgf^ the horizontal Reynolds stress, and a vertical eddy viscosity, 
while the last three terms correspond to the dissipative momentum flux from 
waves to the mean flow, through whitecapping, wave-turbulence interactions, 
and bottom friction. Although the momentum lost by the waves via bot- 
tom friction was shown to eventually end up in the bottom (Longuet-Higgins 
2005), the intermediate acceleration of the mean flow, also known as Eulerian 
streaming, is important for sediment transport, and should be included with a 
vertical profile of T^^'''^ concentrated near the bottom, provided that the wave 
boundary layer is actually resolved in the 3D model (e.g. Walstra et al. 2001). 

The GLM mass conservation writes 

^.« + ^«^0. (44) 

at OXa oz 



where the Jacobian J is the determinant of the coordinate transform matrix 
{6ij + d^i/dxj) from Cartesian coordinates to GLM. (Andrews and Mclntyre 
1978a, eq. (4.2)-(4.4) with p« = p^). 
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2.2 glm2-RANS equations in z-coordinates 



Equations (42) and (44) hold from z — —h to z — ( , which covers the entire 
'GLM water column'. All terms in (42) are defined as GLM averages, except 
for the hydrostatic pressure which does correspond to the Eulerian mean 
position. 

For practical numerical modelling, it is however preferable that the height of 
the water column does not change with the local wave height. We will thus 
transform eq. (42), except for , by correcting for the GLM-induced vertical 
displacements. This will naturally remove the divergence of the GLM flow 
related to J 7^ 1. The GLM vertical displacement ^3 is a generalization of eq. 
(23) 



^^{x,z,t) = j E{k)m 



'^sinh [2k{z + h)] _^ ^sinh^ [k{z + h)] k dua 



2smh\kD) 



sinh (kD) a dz 



dk. 



(45) 



and the Jacobian is J = 1 + J2 + 0{ef). Because the GLM does not induce 
horizontal distortions, a vertical distance dz' — Jdz in GLM corresponds to a 
Cartesian distance dz, giving. 



5rf 



dz 

One may note that 



(46) 



/ 



Jd^ = C^ + /i-eJ(o) = L'. 



(47) 



We now implicitly define the vertical coordinate z* with 



(48) 



Any field 4>{xi, X2, z, t) transforms to (p*{xl, X2, z*, t*) with 



90 


d(f)* 


St 90* 


'dt~ 


dt* 


Sz dz* 


(90 


d(p'' 


Sa d(f)* 


dXa 


dxl 


Sz dz* 



(49) 
(50) 
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dcf) _ 1 90* 
dz Sz dz* 



(51) 



with St, Sz and Sq, the partial derivatives of s with respect to t*, z* and x*, 
respectively. The coordinate transform was built to obtain the following iden- 
tity 



s,J^l + [ef) . 



(52) 



Removing the * superscripts from now on, the mass conservation (44) multi- 
plied by Sz may be written as 



dxa dz ' 



(53) 



where the vertical velocity, 



- U^Sa - St 



W 



1 + 0{e) 
dfjdz ' 



(54) 



is the Lagrangian mass flux through horizontal planes. 

Neglecting terms of order ef and higher, the product of (42) and s^J is re- 
written as. 



d 



+ S' 



shear 



-rs-^ H Aq,, 

OZ 



(55) 



with 



w — J 



- UaSa - St 



^W-Ps + 0ia6l62/k), 



(56) 



the quasi-Eulerian advcction velocity through horizontal planes. From now on 
we shall use exclusively these glm2z-RA'NS equations in z coordinate, with a 
non-divergent GLM velocity field u^. 



Using eq. (25), we may re- write (53) as 



dua did 

— - H = 0. 

dxn, dz 



(57) 
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2.2.1 Surface boundary conditions 



Taking an impermeable boundary, the kinematic boundary condition is given 
by Andrews and Mclntyre (1978a, section 4.2), 

f at .=r. (58) 



It is transformed to z coordinates as 



When the presence of air is considered, it should be noted that the GLM posi- 
tion is discontinuous in the absence of viscosity, because the Stokes corrections 
for C have opposite signs in the air and in the water. This discontinuity arises 
from the discontinuity of the horizontal displacement (air and water wave- 
induced motions are out of phase). A proper treatment would therefore require 
to resolve the viscous boundary layer at the free surface. This question is left 
for further investigation. However, we note that due to the large wind veloc- 
ities and possibly large surface currents unrelated to wave motions, a good 
approximation is given by neglecting the Stokes corrections for the horizontal 
air momentum, 

K = ^a+Pa, (60) 

where the — and -|- exponents refer to the limits when approaching the bound- 
ary from below and above, respectively. 

For the mean horizontal stress, we use the results of Xu and Bowen (1994), 

Ta = Snnna + Snsn3 Sit Z^C (61) 



with S the stress tensor, with normal and shear Sns stresses on the surface, 
generally defined by 

(Oil ' du ■ \ 



with V the kinematic viscosity, and the local unit vector normal to the surface, 
to first order in £i. 
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Taking the Lagrangian mean of (61), one obtains, 

8v r)P _ 

= = C + p^i.^ + p^i.^ at z^C: (64) 



where is the total air-sea momentum flux (the wind stress), as can be 
measured above the wave-perturbed layer (e.g. Drennan et al. 1999). is 
the a component of the wave-supported stress due to surface-slope pressure 
correlations. 



r:=P§- (65) 



■'a 



The second viscous term p^udPa/dz was estimated using the GLM average of 
wave orbital shears (Ardhuin and Jenkins 2006), it is the well-known virtual 
wave stress (e.g. Xu and Bowen 1994, eq. 18). That stress corresponds to 
wave momentum lost due to viscous dissipation, and it can be absorbed into 
the boundary conditions because it is concentrated within a few millimeters 
from the surface (Banner et Peirson 1998). At the base of the viscous layer of 
thickness Sg, (64) yields, using an eddy viscosity K^, 

^a-^a - Pw^^— = PwKz^- at z = -dg. (66) 



2.2.2 Bottom boundary conditions 

The same approach applies to the bottom boundary conditions. The kinematic 
boundary condition writes 

-^ + {ua + Pa)-g^ ^(w + Pa) at z^-h"". (67) 



If an adherence condition is specified at the bottom, which shall be used be- 
low, the bottom boundary condition further simplifies as h — h. It may also 
simplify under the condition that the wave amplitude is not correlated with 
the small scale variations of h, which is not generally the case (e.g. Ardhuin 
and Magne 2007). For the dynamic boundary conditions, pressure-slope cor- 
relations give rise to a partial reflection of waves, that may be represented 
by a scattering stress (e.g. Hara and Mei 1987, Ardhuin and Magne 2007). 
This stress modifies the wave pseudo-momentum without any change of wave 
action (see also Ardhuin 2006). 

The effect of bottom friction is of considerable interest for sediment dynamics 
and deserves special attention. For the sake of simplicity, we shall here use the 
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conduction solution of Longuct-Higgins for a constant viscosity over a flat sea 
bed as given in the appendix to the proceedings of Russel and Osorio (1958). 
We shall briefly consider waves propagating along the x-axis, and we assume 
that the mean current in the wave bottom boundary layer (WBBL) is at most 
of the order of the wave orbital velocity outside of the WBBL. Instead of 
(11)-(16) the orbital wave velocity and displacements near the bottom take 
the form, 



Ui = Uq 



COS ip — e ^ cos{ip — z) 



(6J 



w — 

6 = 

^3 = 



2 

2a; 



2zsmip — sm{ip — z)e ^ + sinip + cos{ip — z)e ^ — cosip (69) 



sin'^ — sin('0 — z)e 



(70) 



2zcosip — cos('0 — z)e ^ + cosip + sin('0 — z)e ^ — sinip (71) 



1 /2 

where ip = kx — ujt is the wave phase, 5f = {2v/uj) ' is the depth scale for 
the boundary layer, z = {z + h)/6f is a non-dimensional vertical coordinate, 
Uq = aa/ sinh(/cD) is the orbital velocity amplitude outside the boundary 
layer. 

Based on these velocities and displacements, the wave pseudo-momentum P, 
is 



Pi = - Cs.iw = ^ 1 + e cos(2^) - 2 cos ze' 



(72) 



This is equal to the Stokes drift — Ui^i^i + ^1,3^3 computed by Longuet- 
Higgins. Besides, the rate of wave energy dissipation induced by bottom fric- 
tion is -S'bfric = PwOJUq/2 giving a bottom friction stress T^^^^^dz — kaSuxid {PwO')- 

Generalizing this approach to a turbulent bottom boundary layer (e.g. Longuet- 
Higgins 2005) one may replace the constant viscosity with a depth-varying 
eddy viscosity. If the wave bottom boundary layer (WBBL) is resolved, will 
also include the momentum lost by waves through bottom friction, as given 
by the depth- integral of T^^''"^. One may estimate P from the vertical proflles 
of the wave orbital velocities and w, and the modifled pressure (38) has 
to be corrected for the change in wave orbital velocities in the WBBL. Many 
WBBL models are available for estimating these wave-induced quantities. 

If the bottom boundary layer is not resolved, on may take the lowest model 
level at the top of the wave boundary layer. The bottom stress may then be 
computed from a parameterization of the bottom roughness ^oa' (e.g. Mathisen 
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and Madsen 1996, 1999), which relates the bottom stress 




(73) 



to the current velocity Ua at the lowest model level z, 



— KjUL^q In 



for z-\-h < 5f. 



(74) 



Then the near-bottom velocity Ua should be taken equal to the Eulerian 
streaming velocity ~ 1.5Pa (see e.g. Marin 2004, for turbulent cases with 
rippled beds). Further, in this case the bottom stress should not inchidc 
the depth integral of T^fi^i^ ^jj-^jg latter remark also applies to depth-integrated 
equations. Indeed, r^^ = J_h^^^ T^^^'^d^; is a flux of momentum into the bot- 
tom due to wave bottom friction, r^^ does not participate in the momentum 
balance that gives rise to a sea level set-down and set-up (Longuet-Higgins 



3 Relations between the present theory and known equations 



3. 1 Depth-integrated GLM for a constant density p, 



Using (59) the mass conservation equation in z coordinates (53) classically 
gives (e.g. PhiUips 1977) 



which is exactly the classic shallow- water mass conservation for constant den- 
sity. 



2005). 




(75) 



dD 
'dt 



(76) 
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with the depth-integrated volume flux vectoiO] M deflned by 



M 



u^dz. 



-h 



(77) 



In the momentum equation, the advection terms may be transformed in flux 
form using mass conservation. However, because some of the original GLM 
advection terms are included in the vortex force, the remaining terms do not 
simplify completely. Using (57) one has. 



Pn 



dt 



+ Uf3 



dUr 

dx 



+ w- 



dz 



+ P:- 



dUa 

dz 



d d d 

^ {PwUa) + {PwUpUa) + ^ [Pw {w + P3) Uc 



Un 



dz 



(78) 



Using (59), (67) and (25), and after integration by parts, these advection terms 
integrate to 



dt dx 



pwUaUpdz 

J 



C 



dMJ^ duAa,.^ fpdu 



V(3 



/3 



^73 



where the zeroth order wave advection velocity is deflned by, 
C 

UAaM^ = J UaPpdz, (80) 



which is equal, at lowest order, to the second term in (8). The wave-induced 
mass transport is the depth-integrated pseudo-momentum, 

C 

M"" = f Pdz. (81) 



-h 



Finally, the quasi-Eulerian volume flux is deflned by M"' = M — M*". 

For terms uniform over the depth [dp^ /dxa and dS'^/dxa) the integral is 
simply the integrand times the depth. 

3 Phillips (1977) uses the notation M instead of M, and M instead of M'". 
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It should be noted that the depth-integrated vortex force involves the advec- 
tion velocity u^, 

/ ea3/3 (/3 + ous) Ppdz = e^sp (/a + ^^a) MJ, (82) 

-h 

with 

^^a = eaa/3 {duAp/dxa - duAa/dxp) . (83) 



The vertical integration of (55) thus yields 



dM^ d 



dt ^ dx 



'/3 




+ eaa/s/sM^ + [p^gC + Pa) 



= -e.s, (/a + n,) M, - UAa^ - -^M, + J P,^d. 

—h 

—h —h 

The source of momentum X^°* is simply the sum of the mean momentum 
fluxes at the top and bottom, and the source of momentum due to diabatic 
wave-mean flow interactions (i.e. breaking and wave-turbulence interactions). 

These equations are very similar to those of Smith (2006, eq. 2.29), our term 
5''' is simply termed J in Smith (2006), and corresponds to Smith's kiD^ . 
The only differences are due to the vertical shear in the current. The advection 
velocity UAa replaces Smith's mean flow velocity. Since UAa is the proper lowest 
order advection velocity for the wave action (Andrews and Mclntyre 1978b), 
this is a simple extension of Smith's result to depth-varying currents. The 
term involving 5'^'^'^^'' is also obviously absent from Smith's equations. The last 
differences in (84) are the last two terms on the second line, but they also 
cancel for a depth-uniform current u^- 



3.2 Equations of McWilliams et al. (2004) 



The approach of MRL04 is in the line of perturbation theories presented by 
Mei (1989) for Eulerian variables and monochromatic waves. Although the 
result of MRL04 corresponds to a particular choice of the relative ordering of 
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small parameters, it is given to a high enough order so that it does cover most 
situations at a lower order. In particular MRL04 have pushed the expansion 
to order ef for some terms because they assumed a ratio a/f^ of order ef, with 
£i the wave slope. This ratio, in practice, may only be attained for relatively 
steep wind waves (developed wind seas and swells generally have slopes of the 
order of 0.05). They also assumed that ~ £2 (the wave envelope varies on 
a scale relatively larger than the wavelength compared to the present theory 
in which Si ~ £2 is possible). These authors also separated the motion into 
waves, long waves and mean flow, and considered in detail the rotational part 
of the wave motion caused by the vertical shear of the current. 

MRL04 thus obtained Eulerian-mean equations that only correspond to mea- 
surable Eulcrian averages under the level of the wave troughs. Because they 
use an analytic continuation of the velocity profiles across the air-sea inter- 
face, the physical interpretation of their average is unclear between the crests 
and troughs of the waves. We shall neglect here their terms of order ef (i.e. 
terms that involve the wave amplitude to the power of four), which amounts 
to choosing a slightly different scaling. Since we shall consider here random 
waves, this avoids cumbersome considerations of the wave bispectrum. 

The Eulerian-mean variables of MRL04 should be related to the Lagrangian 

mean values by the Stokes corrections (5), so that their horizontal Eulerian- 
mean velocity q corresponds to — u'^. Because they have subtracted the 
hydrostatic pressure with the mean water density p^Q, their mean pressure {p) 
should be equal to the Eulerian mean pressure p + Pwogz, with p related to 
the GLM pressure via eq. (37). 

Absorbing the long waves in the mean flow (i.e. allowing the mean flow to vary 
on a the wave group scale, see also Ardhuin et al. 2004), MRL04 equations for 
the 'Eulerian' mean velocity {qi, q2,w) can be written as 



= (Pw - Pwo)9 ~ ^ + ^2) + K 



d 



dqp ^ dw 
dxjs dz 



{p) = P^g[C-kEFscFss)-Vo at z^O 
w = -w"^* at z = 



dXr 



(/Ci+/C2) + J„ 

(85) 
(86) 

(87) 



(89) 



with 
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J, 



K 



a 




(90) 



(91) 



(92) 



/C2 



(93) 



-h 




(94) 



The original notations of MRL04 (sec also Lane et al. 2007) have been trans- 
lated to the notations used above and order ef terms have been neglected. 

These equations are clearly analogue to the glm2z-RANS equations presented 
here. In particular the vertical vortex force term K corresponds to our Ki 
that gets into S'^^'^'^^, the dynamically relevant kinematic pressure pressure 
(p) + /Ci + /C2 corresponds to our pressure p defined by (38), and the verti- 
cal Stokes velocity W^* corresponds to our P3. There are only two differences. 
One is between the surface boundary conditions for these two pressures, with 
a difference only due to }C2{z — 0) — -f^2(C )■ Integrating by parts to es- 
timate ]C2{z — 0), this difference is found to be of the order of gkEs^. Such 
a difference is of the same order as extra terms that would arise when using 
wave kinematics to first order in the current curvature (Kirby and Chen 1989), 
and properly transforming u in u. The second difference between MRL04 and 
the present equations is that the wave pseudo-momentum P differs from the 
Stokes drift u*^ when the current shear is large, and both generally differ from 
the expression for u'^ given by MRL04. Since MRL04 took the current and 
wave orbital velocity to be of the same order, in that context the difference 
P — u'^ is of higher order and thus the two sets of equations are consistent in 
their common range of validity. 

A general comparison of 2D depth-integrated equations is discussed by Lane 
et al. (2006). The present work therefore brings a further verification of their 
3D form of the equations, and an extension to relatively strong currents, pos- 
sibly as large as the phase velocities. As expected, the Eulerian averages of 
McWilliams et al. (2004) arc identical to the quasi- Eulerian fields in GLM 
theory, because they obey the same equations, except for current profile cur- 
vature effects, which were partly neglected here. The "Eulerian" mean current 
of MRL04 can thus be physically interpreted as a quasi-Eulerian average, de- 
fined as the GLM average minus the wave pseudo-momentum. Except for a 
Jacobian that introduces relative corrections of second order in the wave slope, 
this averaging is identical to the procedure used by Swan et al. (2001). Above 
the trough level, this average should not be confused with a truly Eulerian 
average, as obtained from in-situ measurements for example. In such mea- 



26 



surcmcnts the Stokes drift would be recorded in the trough-to-crest region 
(figure l.a). 



4 Limitations of the approximations 

The (7/m2z-RANS equations have been obtained from the exact GLM equa- 
tions, under 6 restricting hypotheses related to the wave slope and Ursell 
number (HI and H2), the horizontal scales of variation of the wave amplitude 
(H3), the current profile (H4 and H5) and the vertical mean velocity (H6). 
These hypotheses essentially allowed us to use the linear wave-induced quan- 
tities given by cqs. (11)-(19). In practical conditions, these hypotheses may 
not be verified and the resulting glm2z-IiANS equations may have to be mod- 
ified. Here we investigate the importance of H3, H2 and HI, using numerical 
solutions from an accurate coupled mode model for irrotational wave propa- 
gation over any bottom topography, and an accurate analytical solution for 
incipient breaking waves, respectively. 

4.1 Bottom slope and standing waves 

In absence of dissipation and given proper lateral boundary conditions the 
flow in wave shoaling over a bottom slope is irrotational and can thus be 
obtained by a numerical exact solution of Laplace's equation with bottom, 
surface, and lateral boundary conditions. For waves of small amplitudes this 
can be provided by a solution to this system of equations to second order in 
the wave slope. Belibassakis and Athanassoulis (2002) have developed a second 
order version of the National Technical University of Athens numerical model 
(NTUA-nl2) to solve this problem in two dimensions. Here we apply their 
model to the simple case of monochromatic, unidirectional waves propagating 
along the x axis, with a topography uniform along the y axis. The topogra- 
phy h{x) varies only for < a; < L and is constant h{x) = /ii for x < 
and h{x) = h2 for a; > L. In that case the Eulerian mean current V0o(x) is 
irrotational, and uniform over the vertical as x approaches ±00 (e.g. Belibas- 
sakis and Athanassoulis 2002, table 1 and flgure 5). We shall further restrict 
our investigation to the case of a monochromatic wave train of known radian 
frequency lu and incident amplitude a, giving rise to reflected and transmitted 
wave trains of amplitudes Ra and Ta. Numerical calculations are given for a 
bottom proflle as given by Roseau (1976) for which the reflection coefficient it! 
is known analytically, thus providing a check on the quality of the numerical 
solution. 

The bottom is deflned here by x and z coordinates given by the real and 
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imaginary part of the complex parametric function of the real variable x\ 

+ + + (96) 



We choose /ii = 6 m and /i2 = 4 m and a wave frequency of 0.19 Hz {uj — 
1.2 rad s~^). For cto = IStt/ISO the maximum bottom slope is £2 = 2.6 x 10~^ 
(figure 1), and the reflection coefficient for wave amplitude is i? = 1.4 x 10~^ 
(Roseau 1976), so that reflected waves may be neglected in the momentum 
balance. Due to the shoaling of the incident waves, the mass transport induced 
by the waves increases in shallow water, and thus the mean current must 
change in the x direction to compensate for the divergence in the wave-induced 
mass transport. We shall further take a zero- mean surface elevation as x — > 
—00. The second order mean elevation is obtained as a result of the model. 
We also verified that the vertical wave pseudo-momentum compensates for the 
divergence of the horizontal component so that in this case for linear waves 
the wave pseudo- momentum is non-divergent (figure 3). 

For mild bottom slopes, the refiection coefficient is small as predicted by 
Roseau (1976). The NTUA-nl2 model used here generally gives accurate refiec- 
tion coefficients, but it tends to overestimate very weak reffections. In the first 
case investigated here, the numerical refiection is i? = 1 x 10~^, with no sig- 
nificant effect on the wave dynamics. The NTUA-nl2 model is used to provide 
the Fourier amplitudes of the mean, ffrst and second harmonic components of 
the velocity potential, over a grid of 401 (horizontal) by 101 (vertical) points. 
From these discretized potential fields, the mean, first and second harmonic 
velocity components are obtained using second order centered finite differ- 
ences. As expected, the numerical solution gives a horizontal mean flow u that 
compensates the divergence of the wave mass transport and is thus of order 
a/ke^. Further u is almost uniform over the vertical and is irrotational (fig- 
ure 2.b). The vertical mean velocity is of higher order. The GLM momentum 
balance is thus dominated by the hydrostatic and dynamic pressure terms 
and S"^ . Although these two terms are individually of the order of 0.01 m^ s~^, 
their sum is less than 2 x 10~^^ m^ s~^ in the entire domain, at the roundoff 
error level. It thus appears that this part of the momentum balance is much 
more accurate than expected from the asymptotic expansion. Indeed, for any 
bottom slope, in the limit of small surface slopes and for irrotational ffow and 
periodic waves, the Stokes correction (5) for the pressure and the time aver- 
age of the Bernoulli equation give the following expression for the modified 
kinematic pressure (38) 



UjUj p ^ 1 ^ dp UjUj 
Pw 2 ^ dxj 2 
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Fig. 2. (a) Instantaneous pressure perturbation (p—p) / (pwO) given by the NTUA-nl2 
model (Belibassakis and Athanassoulis 2002), including the second order Stokes 
component in waves with amplitude a = 0.12 m, over the bottom given by eq. (95). 
(b) Mean current —u, and (c) horizontal wave pseudo-momentum Pi estimated 
from eq. (7), and verified to be equal to the Stokes drift. Arrows indicate the flow 
directions. 
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(96) 
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Fig. 3. Vertical wave pseudo-momentum for the same case as figure 2, estimated 
from eq. (7), and verified to satisfy (25). 

where the equahties only hold to second order in the surface slope. Thus the 
kinematic modified pressure p has no dynamical effect to second order in the 
wave slope, as already discussed by McWilliams et al. (2004) and Lane et al. 
(2007). For irrotational flow, this remains true for any bottom topography 
and even for rapidly varying wave amplitudes, including variations on scales 
shorter than the wavelength. 

Thus the only wave effect is the static change in mean water level (set-up or 
set-down), and dynamic consequences in the WBBL, where S'^ goes to zero, 
leaving the hydrostatic pressure gradient to drive a mean flow that can only 
be balanced by bottom friction. For slowly varying wave amplitudes the mean 
sea level is given by Longuet-Higgins (1967, eq. Fl) 



kE 



+ 



koEo 



smh{2kD) sinh(2A;o 



(97) 



where the subscript correspond to quantities evaluated at any fixed horizon- 
tal position, the choice of which being irrelevant to the estimation of horizontal 
gradients of (. 

Equation (97) is well verified by the NTUA-nl2 result for the case considered 
so far (figure 4. a). However, this is no longueur true for rapid variations in 
the wave amplitude a{x), i.e. due to partially standing waves. In that case one 
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should use Longuet-Higgins' eq. D (op. cit.) 



C{x) = - 



UpUjB - -»3 
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(98) 
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with and given by linear wave theory. Eq. (98) is a generalization of 
Miche's (1944a) mean sea level solution under standing waves. Contrary to 
propagating wave groups, for which the mean sea level is depressed under 
large waves, here the depression occurs at the nodes of the standing wave, 
where the horizontal velocities are largest and amplitudes are smallest (figure 
4.c). 

Eq. (98) is well verified in the presence of partially standing waves. To il- 
lustrate this, we have modified the bottom topography, adding a sinusoidal 
bottom perturbation for x > 180 m with an amplitude of 5 cm and a bot- 
tom wavelength half of the local waves' wavelength, which maximizes wave 
reflection (Kreisel 1949). This yields a wave amplitude reflection R = 0.03, for 
cj = 1.2 rad s~^, of the order of observed wave reflections over gently sloping 
beaches (e.g. Elgar et al. 1994). The bottom is shown on figure 4.b. Although 
the standing wave pattern is hardly noticeable in the surface elevation (the 
amplitude modulation is only 6%, figure 4.c), the small pressure modulation 
occur at much smaller scales, so that the associated gradient can overcome 
the large scale gradients of the hydrostatic pressure (figure 4.d). As a result 
small partial stading waves can dominating the momentum balance in the 
WBBL (see Longuet-Higgins 1953, Yu and Mei 2000 for solutions obtained 
with constant viscosity). 

In the presence of such standing waves, and in the absence of strong wave 
dissipation, the hydrostatic pressure on the scale of the standing waves (e.g. 
given by Miche 1944a) drives the fiow in the WBBL towards the nodes of the 
standing wave (Longuet-Higgins 1953), and is balanced by bottom friction. 
This WBBL flow drives an opposite flow above, closing a secondary circula- 
tion cell. This secondary circulation is important for nearshore sediment trans- 
port just outside of the surf zone (Yu and Mei 2000). If these sub-wavelength 
circulations are to be modelled, the present gl'm2z-RANS theory should be 
extended to resolve the momentum balance on the scale of partial standing 
waves. 

This extension is relatively simple as it only introduces additional stand- 
ing wave terms in all quadratic wave-related quantities, arising from phase- 
couplings of the incident and reflected waves. This extension provides a gen- 
erahzation of eq. (98) in the presence of other processes. For example, eq. (39) 
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Fig. 4. (a) Mean sea level obtained with the NTUA-nl2 model (Belibassakis and 
Athanassoulis 2002) and the theory of Longuet-Higgins (1967 eq. Fl: without stand- 
ing waves) using conservation of the wave energy flux along the profile, (b) modified 
bottom profile resulting in a 3% amplitude reflection at a; = 1.2 rad s~^, (c) re- 
sulting mean sea level and normalized local wave amplitude a, (d) mean sea level 
gradient (d). 



now becomes 



kE(k) 
sinh 2k D 
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dk 
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with i?(k) the amplitude reflection coefficient and 2tp'{'k) is the phase of the 
partial standing waves defined by V-?/'' = k and dtp' /dt = —k-UAt such that 
it is zero at the crest of the incident waves. Note that the integral is over the 
incident wave numbers only (e.g. for wave propagation directions from to 
tt). Similar expressions are easily derived for the other wave forcing terms. 
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4-2 Effects of wave non-linearity 



Deep or intermediate water waves do not break very often in most conditions 
(e.g. Banner et al. 2000, Babanin et al. 2001), thus the particular kinematics of 
breaking or very steep waves hkely contributes httle to the average forcing of 
the current. However, most of the waves break in the surf zone and deviations 
from Airy wave kinematics may introduce a systematic bias when the glm2z- 
RANS equations are apphed in that context. Many wave theories have been 
developed that are generally more accurate than the Airy wave theory (e.g. 
Dean 1970). However, they may lack some realistic features found in breaking 
waves, such as sharp crests. In order to explore the magnitude of this bias, we 
shall use the kinematics of two-dimensional incipient breaking waves as given 
by the approximate theory of Miche (1944b). 

Miche's theory is based on the asymptotic expansion of the potential flow 
from the triangular crest of a steady breaking wave, extending Stokes' 120° 
corner flow to finite depth. From this Miche obtained his criterion for the max- 
imum steepness of a steady breaking wave, i.e. h/X = 0.14tanh(A;/i) with h the 
breaking wave height and A the wavelength, which favorably compares with 
observations. The Miche wave potential and streamfunction ip are expressed 
implicitly as a function G of the coordinates x — Xc + i{z — Zc), with origin on 
the wave crest {xc, z^)- The coefficients in the series representing the reciprocal 
function G' are obtained from the boundary condition at the surface and bot- 
tom. Unfortunately, these are imposed only under the wave crest and trough, 
so that the bottom streamline may not be horizontal away from the crest. 
This is particularly true for small values of kh. Due to the expansion of G' in 
powers of + 1?/^, the shape of the wave is nevertheless accurate near the crest, 
and since the overall drift velocities are dominated by the corner flow near the 
crest (see also Longuet-Higgins 1979), the approximations of Miche have little 
consequence on the drift velocities. The function G' was modified here to make 
the bottom actually fiat, and the vertical under the trough an equipotential. 
This deformation adds a weak rotational component to the motion and the 
wave streamlines are weakly modified at the bottom under the wave troug 10. 
The resulting wave for kh = 0.58 (corresponding to 6 = 1 in Miche 1944b) 
is shown in figure 5. a. A numerical evaluation of that solution is obtained at 
201 equally spaced values of ip and 401 equally spaced values of (figure 5.b). 
The GLM displacement field ^ is computed as described in section 2.1. Since 
the streamlines are known in the frame of reference of the wave, Lagrangian 
positions of 201 particles initially placed below the crest at Xi{0) = 0, were 
tracked over four Eulerian wave periods. The positions {xi{t) , Zi{t)) are given 
by the potential 0j(t) and streamfunction ipi. The Lagrangian period for each 



^ This correction leads to negligible differences compared to the exact solution as 
verified with streamfunction theory to 60th order. 
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particle T^" is determined by detecting the first time when the particles pass 
under the crest again. The Lagrangian mean velocity of each particle is then 

Xi{T^)/T{", and it corresponds to a vertical position Zi = /g ' Zi{t)dt. This 
defines the Lagrangian mean velocity u^{zi) in GLM coordinates. Following 
the coordinate transformation in section 2, we further transform the GLM 
velocity profile to z coordinate (figure 5.c). The resulting profile of u-^ has a 
horizontal tangent at z = 0, as discussed by Miche (1944b). 

Contrary to Miche (1944b) who defined the phase speed C of his wave by 
imposing a zero mass transport, we have defined C so that P = u"^ with the 
pseudo-momentum P estimated from eq. (7) using finite differences applied to 
the displacement field. The two profiles of P, estimated from eq. (7), and U''", 
estimated by time integration of particle positions coincide almost perfectly. 
Thus the estimation of P provides a practical method for separating the mean 
current from the wave motion. Starting from any value of C, the difference 
between and P is the mean current velocity u. Here C was corrected to 
have u = 0. 

Prom ^, Bernoulli's equation can be used to obtain the GLM of velocities 
and pressure. Compared to linear wave theory, the Stokes drift in a Miche 

wave is much more sheared. It should be noted that in the cnoidal theory in- 
vestigated by Wicgel (1959) this drift velocity is depth- uniform. Thus cnoidal 
wave theories may produce inaccurate results for 3D wave-current interactions 
when extrapolated to breaking waves. This marked difference in the 3D mean 
fiow forcing due to breaking waves compared to linear waves calls for a deeper 
investigation of this question. Investigating such kinematics, may provide a ra- 
tionale for the parameterization of nonlinearity in the glm2z-'RANS equations 
proposed here. Such a parameterization is proposed by Rascle and Ardhuin 
(manuscript in preparation for the Journal of Geophysical Research). 



5 Conclusion 

We have approximated the exact Generahzed Lagrangian Mean (GLM) wave- 
averaged momentum equations of Andrews and Mclntyre (1978a), to second 
order in the wave slope, allowing for strong and sheared mean currents with 
limited curvature in the current profile. These approximated equations were 
then transformed by a change of the vertical coordinate, giving a non-divergent 
GLM fiow in z coordinates. The resulting conservation equations for horizontal 
momentum (55) and mass (57), with boundary conditions (59)-(74) may be 
solved using slightly modified versions of existing primitive equations models, 
forced with the results of spectral wave models. Although the Stokes drift 
introduces a source of mass at the surface for the quasi- Eulerian flow, this 
is does not pose any particular problem, and such mass source have long 
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Fig. 5. (a) Illustration of the drift over 2 Eulerian periods in periodic Miche waves. 
Trajectories are color-coded with their initial depth, below a wave crest. The thin 
black lines are the lines of constant potential and streamfunction at t = 0. (b) Field 
of displacements defining the GLM, as in figure I.e. The dash-dotted line is the 
GLM position of the free surface ( . (c) Profiles of Eulerian and Lagrangian mass 
transport velocity in a Miche wave compared to a linear wave with the same values 
of k and h. 



been introduced for the simulation of upwellings. The HYCOM model (Bleck 
2002) was modified by R. Baraille to solve a simplified set of the present 
equations, retaining only the wave-induced mass transport in both the mass 
and momentum equations, and the tracer equation (in which the advection 
velocity is simply u^, see also MRL04). This work was applied to the a hindcast 
of the trajectories of sub-surface oil pellets released by the tanker Prestige- 
Nassau, which sank off Northwest Spain in November 2002 (presentation at 
the 2004 WMO-JCOMM 'Oceanops' conference held in Toulouse, France). 
The full equations derived here have also been implemented in the ocean 
circulation model ROMS (Shchepetkin and Mc Williams 2003), and results will 
be reported elsewhere. The equations presented here have also been applied 
for the modelling of the ocean mixed layer in horizontally- uniform conditions 
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(Rascle et al. 2006). 



Although a general expression for the turbulent closure has been given, it has 
not been made explicit in terms of the wave and mean flow quantities beyond a 
heuristic closure that combines an eddy viscosity mixing term with the known 
sources of momentum due to wave dissipation. A proper turbulent closure 
is left for further work, possibly extending and combining the approaches of 
Groeneweg and Klopman (1998), with those of Teixeira and Belcher (2002). 
Further, some wave forcing quantities have been expressed in terms of the 
Eulerian mean current u instead of the quasi-Eulerian mean current u. The 
conversion from one to the other, can be done using eq. (24), to the order of 
approximation used here. However, it would be more appropriate, in particular 
for large current shears, to start from quasi-Eulerian wave kinematics, instead 
of Eulerian solutions of the kind given by Kirby and Chen (1989, our eq. 
10-12). 

Beyond the turbulence closure, there are essentially two practical limitations 
to the approximate glm2z-IiANS equations derived here. First, the expansion 
of wave quantities to second order in the surface slope is only qualitative in the 
surf zone. Although this was acceptable in two dimensions (see Bowen 1969 
and most of the literature on this subject), it is expected to be insufficient 
in three dimensions due to a significant difference in the profile of the wave- 
induced drift velocity P, which exhibits a vertical variation with surface values 
exceeding bottom values by a factor of 3, even for kh < 0.2 in which case linear 
wave theory predicts a depth- uniform P. This conclusion is based on both the 
approximate theory of Miche (1944b), and results of the streamfunction theory 
of Dalrymple (1974) to 80th order. Such numerical results can be used to 
provide a parameterization of these effects. Further investigations using more 
realistic depictions of the kinematics of breaking waves will be needed. Second, 
the vertical profile of the mean current in the surf zone may be such that the 
wave kinematics are not well described by the approximations used here. A 
strong nonlinearity combined with a strong current shear and curvature can 
lead to markedly different wave kinematics (e.g. da Silva and Peregrine 1988). 

With these caveats, the equations derived here provide a generalization of 
existing equations, extending Smith (2006) to three dimensions and vertically 
sheared currents, or McWilliams et al. (2004) to strong currents. Of course, 
mean flow equations can be obtained, at least numerically, using any solution 
for the wave kinematics with the original exact GLM equations, as illustrated 
in section 4.2. The wave-forcing on the mean flow is a vortex force plus a 
modified pressure, a decomposition that allows a clearer understanding of the 
wave-current interactions, compared to the more traditional radiation stress 
form. This is most important for the three-dimensional momentum balance 
and/or in the presence of strong currents, e.g. when a rip current is widened 
by opposing waves, as observed by Ismail and Wiegel (1983) in the laboratory. 
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Such a situation was also recently modelled by Shi et al. (2006). 
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Symbol 



name 



where defined 



i and I 


indices of tfie horizontal dimensions 


alter (oj 


q 
o 


index of the vertical dimension 


arter (p) 


a 


wave amplitude 


after (12) 


U = h^-C, 


mean water depth 


alter yt ) 


t - [J1,J2,J3) 


Coriolis parameter vector (twice the rotation vector) 


arter (p ) 


TP TP TP A TP 

i'cc, ^cs, ^sc and l<ss 


Vertical profile functions 


alter (Izj 


9 


acceleration due to gravity and Earth rotation 


alter 1^1 ) 


u 
a 


depth of the bottom (bottom elevation is z = —h) 


before (8) 


J 


Jacobian of GLM average 


alter (^44j 


k = (A;i,fc2j 


wavenumber vector 


alter (7j 


T/' 


Depth-integrated vertical vortex force 


loo\ 

(33) 


K2 


Shear-induced correction to Bernoulli head 


(29) 




vertical eddy viscosity 


(43) 


{■y 


Lagrangian perturbation 


(2) 




Lagrangian mean 


(1) 


m 


shear correction parameter 


(20) 


M 


depth-integrated momentum vector 


(77) 


M"' 


depth-integrated wave pseudo-momentum vector 


(81) 


M"^ 


depth-integrated mean flow momentum vector 


after (81) 


n 


unit normal vector 


(63) 
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Symbol 


name 


where defined 




P 


full dynamic pressure 


after (26) 




P 


wave-induced pressure 


(10) 






hydrostatic pressure 


after (35) 


p = 




wave pseudo-momentum 


(6) 




t 


time 


before (1) 


u = 


■■ {ui,U2,U3) 


velocity vector 






U 


wave-induced velocity 


(11) and (68) 






Lagrangian mean velocity 


after (1) 






advection velocity for the wave action 


(80) 


Ua 


= n^-Pa 


quasi-Eulerian horizontal velocity 


before (24) 


S 


T 


GLM to z transformation function 


(48) 






Stokes correction 


(5) 




Sij 


stress tensor 


(62) 






wave-induced kinematic pressure 


(39) 




^Shear 


shear-induced correction to S'^ 


(40) 




W = U3 


vertical velocity 


before (30) 


W 


= ^3-P3 


quasi-Eulerian vertical velocity 


before (30) 




w 


GLM vertical velocity in z coordinates 


(54) 


X = 


■ {xi,X2,X3) 


position vector 


before (1) 
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Symbol name where defined 

X diabatic source of momentum after (24) 

X diabatic source of quasi-Eulerian mean momentum (27) 

z = X2, vertical position after (8) 

a and /? dummy indices for horizontal dimensions 

(5ij Kronecker's symbol, zero unless i = j after (26) 

£ generic small parameter after (8) 

£\ maximum wave slope after (7) 

£2 maximum horizontal gradient parameter after (7) 

£3 maximum current curvature parameter (9) 

CijkAjBk component i of the vector product A x B after (6) 

free surface elevation before (8) 

A wavelength section 4.2 

u kinematic viscosity of water after (62) 

^ = (^1,^2,^3) wave-induced displacement before (1) 

Pyj density of water (constant) after (12) 

a relative radian frequency after (7) 

Tij mean stress tensor (61) 

ijj wave phase after (7) 

uj absolute radian frequency after (7) and (8) 

^^3 depth-weighted vertical vorticity of the mean flow (83) 

V horizontal gradient operator after (7) 
Table 2 
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